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Abstract The aim of this work is to review and also explore 
even further the escape properties of orbits in a dynamical 
system of a two-dimensional perturbed harmonic oscillator, 
which is a characteristic example of open Hamiltonian sys¬ 
tems. In particular, we conduct a thorough numerical investi¬ 
gation distinguishing between trapped (ordered and chaotic) 
and escaping orbits, considering only unbounded motion for 
several energy levels. It is of particular interest, to locate the 
basins of escape towards the different escape channels and 
connect them with the corresponding escape periods of the 
orbits. We split our examination into three different cases 
depending on the function of the perturbation term which 
determines the number of escape channels on the physical 
space. In every case, we computed extensive samples of or¬ 
bits in both the physical and the phase space by integrat¬ 
ing numerically the equations of motion as well as the vari¬ 
ational equations. In an attempt to determine the regular 
or chaotic nature of trapped motion, we applied the SALI 
method as a chaos detector. It was found, that in all studied 
cases regions of trapped orbits coexist with several basins 
of escape. It was also observed, that for energy levels very 
close to the escape value the escape times of orbits are large, 
while for values of energy much higher than the escape en¬ 
ergy the vast majority of orbits escape very quickly or even 
immediately to infinity. The larger escape periods have been 
measured for orbits with initial conditions in the boundaries 
of the escape basins and also in the vicinity of the fractal 
structure. Most of the current outcomes have been compared 
with previous related work. We hope that our results will be 
useful for a further understanding of the escape mechanism 
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of orbits in open Hamiltonian systems with two degrees of 
freedom. 

Keywords Hamiltonian systems; harmonic oscillators; 
numerical simulations; escapes; fractals 

1 Introduction 

Escaping particles from dynamical systems is a subject to 
which has been devoted many studies over the years. Espe¬ 
cially the issue of escapes in Hamiltonian systems is directly 
related to the problem of chaotic scattering which has been 
an active field of research over the last decades and it still 
remains open (e.g., (8|-[l0|[T9l|20||23l[29l[37l|4ll-|SllSQ|- 
[54]). It is well known, that particular types of Hamiltonian 
systems have a finite energy of escape and for lower values 
the equipotential surfaces of the systems are close and there¬ 
fore escape is impossible. For energy levels beyond the es¬ 
cape energy however, these surfaces open creating exit chan¬ 
nels through which the particles can escape to infinity. The 
literature is replete with studies of such “open" Hamiltonian 
systems (e.g., QE31|2S1BZI0SI|3|-I2I)- 

Usually, the infinity acts as an attractor for an escape 
particle, which may escape through different channels (ex¬ 
its) on the equipotential curve or on the equipotential sur¬ 
face depending whether the dynamical system has two or 
three degrees of freedom, respectively. Therefore, it is quite 
possible to obtain basins of escape, similar to basins of at¬ 
traction in dissipative systems or even the Newton-Raphson 
fractal structures. Basins of escape have been studied in sev¬ 
eral papers (e.g., [ ElEIllffl®])- The reader can find more 
details regarding basins of escape in |2H . The key idea of 
studying escape of orbits in open dynamical systems is the 
existence of a chaotic invariant set of orbits embedded in the 
system and its stable and unstable manifold, where the un¬ 
stable manifold in fact coincides with the fractal boundary. 
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One of the most characteristic models for time-independent 
Hamiltonian systems of two degrees of freedom is undoubt¬ 
edly the well-known Henon-Heiles system [ 36 ). A huge load 
of research on the escape properties of this system has been 
conducted over the years (e.g., [Q]- 3|[6l[]j][^[49]). At this 
point we should emphasize, that all the above-mentioned 
references on escapes in the Henon-Heiles system are ex¬ 
emplary rather than exhaustive, taking into account that a 
vast quantity of related literature exists. 

During the last half century, dynamical systems made 
up of perturbed harmonic oscillators have been extensively 
used in order to describe local motion (i.e., near an equi¬ 
librium point) (e.g., J5l H21 [14] - H6l [33] - [361 [471 EH [6T]- 
[63)). In an attempt to reveal and understand the nature of 
orbits in these systems, scientists have used either numeri¬ 
cal (e.g., n~6l 39 64|) or analytical methods (e.g., fT3li27. 
ISDESBED). Furthermore, potentials made up of harmonic 
oscillators are frequently used in galactic Astronomy, as a 
first step for distinguishing between ordered and chaotic lo¬ 
cal motion in galaxies, since it is widely accepted that the 
motion of stars near the central region of a galaxy can be 
approximated by harmonic oscillations. 

A simple dynamical system of two coupled harmonic 
oscillators for various values of the energy above the escape 
energy has been investigated in (22), where it was found that 
stable periodic orbits are surrounded by stability islands that 
never escape. A further numerical analysis of the same dy¬ 
namical system in ED, revealed that as the energy increases 
beyond the escape value, the majority of chaotic orbits es¬ 
cape either directly, or after a small or large number of in¬ 
tersections with the y = 0 axis. In the same vein, the effects 
of different types of perturbations on both the topology and 
the escaping dynamics in the Henon-Heiles system was ex¬ 
amined in m , where basins of escape were found to exist 
in the physical (x,y) as well as in the phase (y,y) space. 

Escaping and trapped orbits in stellar systems is an is¬ 
sue of paramount importance. In a recent article (60), we 
explored the nature of the orbits of stars in a galactic-type 
potential, which can be considered to describe local motion 
in the meridional plane (R,z) near the central parts of an ax¬ 
ially symmetric galaxy. It was observed, that apart from the 
trapped orbits there are two types of escaping orbits, those 
which escape fast and those which need to spend vast time 
intervals inside the limiting curve before they find the exit 
and eventually escape. Furthermore, the chaotic dynamics 
within a star cluster embedded in the tidal field of a galaxy 
was explored in E). In particular, by scanning thoroughly 
the phase space and obtaining the basins of escape with the 
respective escape times it was revealed, that the higher es¬ 
cape times correspond to initial conditions of orbits near the 
fractal basin boundaries. 

Thus, taking into account all the above-mentioned facts, 
we decided to use a potential of a perturbed harmonic os¬ 


cillator with such perturbing terms producing between two 
and four escape channels in the physical (x,y) space. Here 
we must point out, that these dynamical systems have been 
studied thoroughly in many previous papers so, in the cur¬ 
rent work, we shall try to review the main properties of them 
and also present some more detailed results regarding the es¬ 
cape mechanism of orbits. The aim of this work, is twofold: 

(i) to distinguish between trapped and escaping orbits and 

(ii) to locate the basins of escape leading to different escape 
channels and try to connect them with the corresponding es¬ 
cape times of the orbits. In the forthcoming Part II, we will 
consider open Hamiltonian systems with n(n> 5) channels 
of escape which however, have not been explored yet. 

The present article is organized as follows: in Section [2] 
we describe the properties of the potential we chose for our 
investigation of trapped and escaping orbits. The computa¬ 
tional methods used in order to determine the nature (or¬ 
dered/chaotic and trapped/escaping) of orbits are described 
in Section [3] In the following Section, we conduct a thor¬ 
ough analysis of several sets of initial conditions of orbits 
presenting in detail all the numerical results of our computa¬ 
tions. Our article ends with Section^ where the conclusions 
and the discussion of this research are presented. 


2 Properties of the model potential 


The general form of a two-dimensional perturbed harmonic 
oscillator is 

V(x,y) = 1 (<o 2 x 2 + ©fy 2 ) +eVi(x,y), (1) 

where C0\ and CO 2 are the unperturbed frequencies of oscilla¬ 
tions along the x and y axes respectively, £ is the perturbation 
parameter, while V\ is the function containing the perturbing 
terms. This is called a two-dimensional perturbed elliptic os¬ 
cillator. 

In the present paper, we shall use a two-dimensional per¬ 
turbed harmonic oscillator at the 1:1 resonance, that is when 
0)i = 0)2 = 0), in order to investigate the escape properties 
of orbits. The corresponding potential is 


0) 2 

V{x,y) = — (x z +y 2 ) +eVi(x,y), (2) 

being 0) the common frequency of oscillations along the two 
axes. Without the loss of generality, we may set 0) = 1 and 
£ = 1 for more convenient numerical computations. 

The basic equations of motion for a test particle with a 
unit mass are 



dV 

~dy’ 


(3) 


where, as usual, the dot indicates derivative with respect to 
the time. Furthermore, the variational equations governing 
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the evolution of a deviation vectoiQw = (ii.v. Sy, Sx, Sy) are 

(Sx) = Sx, (Sy) = Sy, 


d 2 V 

dxdy 


Sy, 


d 1 2 V 

(Sx) = 

ax L 

d 2 V d 2 V 


(4) 


The Hamiltonian to potential ^ (with CO = £ = 1) reads 

H=^(x 2 * +y 2 +x 2 +y 2 )+V l (x, y)=h, (5) 


where x and y are the momenta per unit mass conjugate to v 
and y respectively, while h > 0 is the numerical value of the 
Hamiltonian, which is conserved. The Hamiltonian can also 
be written in the form 


H = H 0 +H U 


( 6 ) 


with Hq being the integrable term and H\ the non-integrable 
correction. 

Potential |2]) has a finite energy of escape (h esc ) which 
can be derived as follows: First we solve the system 


dV „ dV „ 
!h~ 0, ~dy~ 0 ' 


(7) 


The solutions of system Q gives all the critical points of 
potential function. The saddle points of ^ are those of the 
critical points that satisfy the condition 


n fd 2 V\fd 2 V\ ( d 2 V \ ^ 

s ~ \~dx 2 ) \~dy 2 )~\d^J <0 ' 8 

The value of the escape energy is obtained, if we insert the 
solution of system 0 which satisfy the condition 0 in 
the potential 0. It becomes evident, that the escape energy 
strongly depends on the particular function of the perturba¬ 
tion term V\ (v, y ). Here we should note, that in the case when 
more than one solutions satisfy simultaneously the condition 
0, then h esc is the minimum of the corresponding values of 
h that are calculated. 


3 Computational methods 

In order to study the escape process in our Hamiltonian sys¬ 
tem, we need to define samples of orbits whose properties 
(escaping or trapped) will be identified. The best method for 
this purpose, would have been to choose the sets of initial 
conditions of the orbits from a distribution function of the 
system. This however, is not available so, we define for each 
set of values of the energy (all tested energy levels are above 
the escape energy), dense grids of initial conditions regu¬ 
larly distributed in the area allowed by the value of the en¬ 
ergy. Our investigation takes place both in the physical (x,y) 

1 If S is the 2N dimensional phase space where the orbits of a dy¬ 

namical system evolve on, a deviation vector w, which describes a 

small perturbation of a specific orbit x, evolves on a 2N dimensional 

space T X S tangent to S. 


and the phase (x,x) space for a better understanding of the 
escape mechanism. In both cases, the step separation of the 
initial conditions along the x and y and x and x axes (in other 
words the density of the grid) was controlled in such a way 
that always there are about 50000 orbits (maximum a grid 
of 225 x 225 equally spaced initial conditions of orbits). 

For each initial condition, we integrated the equations 
of motion 0 as well as the variational equations 0 us¬ 
ing a double precision Bulirsch-Stoer FORTRAN algorithm 
(e.g., m ) with a small time step of order of 10 -2 , which 
is sufficient enough for the desired accuracy of our compu¬ 
tations (i.e. our results practically do not change by halv¬ 
ing the time step). Our previous experience suggests that the 
Bulirsch-Stoer integrator is both faster and more accurate 
than a double precision Runge-Kutta-Fehlberg algorithm of 
order 7 with Cash-Karp coefficients. In all cases, the en¬ 
ergy integral (Eq. ([5])) was conserved better than one part 
in 10“ 10 , although for most orbits it was better than one part 
in KT 11 . 

An issue of paramount importance is the determination 
of the position as well as the time at which an orbit es¬ 
capes. When the value of the energy h is smaller than the 
escape energy, the Zero Velocity Curves (ZVCs) are closed. 
On the other hand, when h > h esc the equipotential curves 
are open and extend to infinity. An open ZVC consists of 
several branches forming channels through which an orbit 
can escape to infinity. At every opening there is a highly un¬ 
stable periodic orbit close to the line of maximum potential 
m which is called a Lyapunov orbit. Such an orbit reaches 
the ZVC, on both sides of the opening and returns along 
the same path thus, connecting two opposite branches of the 
ZVC. Lyapunov orbits are very important for the escapes 
from the system, since if an orbit intersects any one of these 
orbits with velocity pointing outwards moves always out¬ 
wards and eventually escapes from the system without any 
further intersections with the surface of section (see e.g., 
l20l ). The passage of orbits through Lyapunov orbits and 
their subsequent escape to infinity is the most conspicuous 
aspect of the transport, but crucial features of the bulk flow, 
especially at late times, appear to be controlled by diffusion 
through cantori, which can trap orbits far vary long time pe¬ 
riods. 

In our computations, we set 10 5 time units as a maxi¬ 
mum time of numerical integration. Our previous experience 
in this subject indicates, that usually orbits need consider¬ 
able less time to find one of the exits in the limiting curve 
and eventually escape from the system (obviously, the nu¬ 
merical integration is effectively ended when an orbit passes 
through one of the escape channels and intersects one of the 
unstable Lyapunov orbits). Nevertheless, we decided to use 
such a vast integration time just to be sure that all orbits have 
enough time in order to escape. Remember, that there are the 
so called “sticky orbits" which behave as regular ones and 
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their true chaotic character is revealed only after long time 
intervals of numerical integration. Here we should clarify, 
that orbits which do not escape after a numerical integration 
of 10 5 time units are considered as non-escaping or trapped. 

The physical and the phase space are divided into the es¬ 
caping and non-escaping (trapped) space. Usually, the vast 
majority of the trapped space is occupied by initial con¬ 
ditions of regular orbits forming stability islands where a 
third integral is present. In many systems however, trapped 
chaotic orbits have also been observed. Therefore, we de¬ 
cided to distinguish between regular and chaotic trapped or¬ 
bits. Over the years, several chaos indicators have been de¬ 
veloped in order to determine the character of orbits. In our 
case, we chose to use the Smaller ALingment Index (SALI) 
method. The SALI l58l has been proved a very fast, reliable 
and effective tool, which is defined as 

SALI(t) =min(d_,d+), (9) 

where d = ||wi(f) — W 2 (*)|| andd + = ||wi(f) + W 2 </)|| are 
the alignments indices, while wi (t) and W 2 (t), are two de¬ 
viations vectors which initially point in two random direc¬ 
tions. For distinguishing between ordered and chaotic mo¬ 
tion, all we have to do is to compute the SALI along time 
interval t max of numerical integration. In particular, we track 
simultaneously the time-evolution of the main orbit itself as 
well as the two deviation vectors wi (t) and Vf 2 {t) in order to 
compute the SALI. The variational equations ([?]), as usual, 
are used for the evolution and computation of the deviation 
vectors. 

The time-evolution of SALI strongly depends on the na¬ 
ture of the computed orbit since when the orbit is regular 
the SALI exhibits small fluctuations around non zero val¬ 
ues, while on the other hand, in the case of chaotic orbits the 
SALI after a small transient period it tends exponentially 
to zero approaching the limit of the accuracy of the com¬ 
puter (10 -16 ). Therefore, the particular time-evolution of 
the SALI allow us to distinguish fast and safely between reg¬ 
ular and chaotic motion (e.g., l65lh . Nevertheless, we have 
to define a specific numerical threshold value for determin¬ 
ing the transition from regularity to chaos. After conducting 
extensive numerical experiments, integrating many sets of 
orbits, we conclude that a safe threshold value for the SALI 
is the value 10 -7 . In order to decide whether an orbit is reg¬ 
ular or chaotic, one may use the usual method according to 
which we check after a certain and predefined time interval 
of numerical integration, if the value of SALI has become 
less than the established threshold value. Therefore, if SALI 
< 10 -7 the orbit is chaotic, while if SALI > 10 -7 the orbit 
is regular. For the computation of SALI we used the LP-VI 
code fl71 . a fully operational code which efficiently com¬ 
putes a suite of many chaos indicators for dynamical sys¬ 
tems in any number of dimensions. 


4 Numerical results 

Our main objective is to distinguish between trapped and es¬ 
caping orbits for values of energy larger than the escape en¬ 
ergy where the Zero Velocity Curves are open and several 
channels of escape are present. Moreover, two additional 
properties of the orbits will be examined: (i) the directions 
or channels through which the particles escape and (ii) the 
time-scale of the escapes (we shall also use the term escape 
period). In the present paper, we explore these aspects for 
various values of the energy h , as well as for three differ¬ 
ent types of perturbation. The function of the perturbation 
term plays a key role as it determines the location as well 
as the number of the escape channels both in the physical 
and the phase space. In particular, three different cases of 
perturbation are considered which produces two, three and 
four channels of escape at the physical (x,y) space respec¬ 
tively. In both cases, the grids of initial conditions of orbits 
whose properties will be examined are defined as follows: 
For the physical (x,y) space we consider orbits with ini¬ 
tial conditions (vo,yo) x o = while the initial value 
of y'o is always obtained from the energy integral 0 as y'o = 
y(xo,xo,h) > 0. Similarly, for the phase (x,x) space we con¬ 
sider orbits with initial conditions (xq,xq) with yo = 0, while 
again the initial value of y'o is obtained from the energy in¬ 
tegral ([5]). 

4.1 Case I: Two channels of escape 

In this case, the perturbation term is Vj (jc,y) = — xy 2 and the 
corresponding Hamiltonian is 

H i = ^ (x 2 + y 2 + * 2 +y 2 ) -xy 2 = h. (10) 

This Hamiltonian system has an escape energy which equals 
to 1/8 and it has been studied extensively in numerous pre¬ 
vious papers (e.g., I20ll22ll25ll38l 55|). This dynamical sys¬ 
tem has a special symmetry; H\ is symmetric with respect 
to y —y —y. The equipotential curves of the potential 0 for 
various values of the energy h are shown in Fig. [1^. The 
equipotential corresponding to the energy of escape h esc is 
plotted with red color in the same plot. The open ZVC at the 
physical (v,y) plane when hm 0.15 > h esc is presented with 
green color in Fig. []J> and the two channels of escape are 
shown. In the same plot, we denote the two unstable Lya¬ 
punov orbits by L\ and L 2 using red color. In Fig. [lJ we 
depict with different colors two orbits, one escaping from 
channel 1 and the other from channel 2, when h = 0.15. 

In our investigation, we shall deal only with unbounded 
motion of test particles for values of energy in the set h = 
{0.13,0.15,0.17,0.19,0.21,0.23,0.25,0.27,0.30}. First of all, 
we will explore the escape process in the physical (x,y) 
plane. Fig. [2] shows the structure of the (x,y) plane for dif¬ 
ferent values of the energy. Each initial condition is colored 
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Fig. 1 (a): Equipotential curves of the potential |2J for various values of the energy h when V\ (x,y) = —xy 2 . The equipotential curve corresponding 
to the energy of escape is shown with red color; (b): The open ZVC at the physical (x,y) plane when h = 0.15. L\ and L 2 indicate the two unstable 
Lyapunov orbits plotted in red; (c): Two escaping orbits when h = 0.15. The orbit which escapes from channel 1 is potted with green color, while 
red color is used for the orbits which escapes through channel 2. 



Fig. 8 Orbital structure of the (jt,ft) -plane when two channels of es¬ 
cape are present. This diagram gives a detailed analysis of the evolu¬ 
tion of the trapped and escaping orbits of the dynamical system when 
the parameter ft changes. The color code is as in Fig. [5] 


according to the escape channel through which the particu¬ 
lar orbit escapes. The gray regions on the other hand, denote 
initial conditions where the test particles do not escape. The 
outermost black solid line is the Zero Velocity Curve (limit¬ 
ing curve) which is defined as V(x,y) = ft. It is seen, that for 
values of energy larger but yet very close to the escape en¬ 
ergy (ft < 0.16) a large portion of the (x,y) plane is covered 
by stability islands which correspond to initial conditions 
of trapped orbits surrounding by a very rich fractal struc¬ 
ture. Looking carefully the grids we also observe that there 
is a highly sensitive dependence of the escape process on the 


initial conditions, that is, a slight change in the initial condi¬ 
tions makes the test particle escape through another channel, 
which is is a classical indication of chaos. As the value of 
the energy increases the stability islands with trapped regu¬ 
lar orbits are reduced and basins of escape emerge. Indeed, 
when ft = 0.30 all the computed orbits of the grid escape and 
there is no indication of bounded motion or whatsoever. By 
the term basin of escape, we refer to a set of initial condi¬ 
tions that corresponds to a certain escape channel. The es¬ 
cape basins become smoother and more well-defined as the 
energy increases and the degree of fractility decrease^] The 
fractility is strongly related with the unpredictability in the 
evolution of a dynamical system. In our case, it can be inter¬ 
preted that for high enough energy levels, the test particles 
escape very fast from the scattering region and therefore, the 
system’s predictability increases. 

Fig. [3^ shows the evolution of the percentages of trapped 
and escaping orbits on the physical (x,y) plane when the 
value of the energy ft varies. One may observes, that when 
ft = 0.13, that is just above the escape energy, trapped, escap¬ 
ing through channel 1 and escaping through channel 2 orbits 
almost share the entire plane. As the value of the energy in¬ 
creases however, the rate of trapped orbits drops rapidly and 
when ft > 0.28 it vanishes. At the same time, the percent¬ 
age of orbits escaping through channel 1 increases steadily 
and for ft > 0.35, it seems to saturate around 65%, thus oc¬ 
cupying around two thirds of the (v,y) plane. On the other 
hand, the rate of orbits escaping through channel 2 increases 
for ft < 0.17 but then is exhibits a slow reduction and for 
ft > 0.37 it saturates around 35%. Therefore, one may con¬ 
cludes that for high energy levels (ft > 0.35), all orbits in 


2 The fat-fractal exponent increases, approaching the value 1 which 
means no fractal geometry, when the energy of the system is high 
enough (see (61). 
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Fig. 2 The structure of the physical (x,y) plane for several values of the energy h, distinguishing between different escape channels. The color 
code is as follows: Trapped (gray); escape through channel 1 (green); escape through channel 2 (red). 


the (x,y) plane escape and about two thirds of them choose 
channel 1. 

The following Fig. [4] shows how the escape times t esc 
of orbits are distributed on the (x,y) plane. Light reddish 
colors correspond to fast escaping orbits, dark blue/purpe 
colors indicate large escape periods, while gray color de¬ 
note trapped orbits. We observe, that when h = 0.13, that is 
a value of energy very close to the escape energy, the es¬ 
cape periods of the majority of orbits are huge correspond¬ 
ing to tens of thousands of time units. This however, is antic¬ 
ipated because in this case the width of the escape channels 
is very small and therefore, the orbits should spend much 
time inside the equipotential curve until they find one of the 
openings and eventually escape to infinity. As the value of 


the energy increases however, the escape channels become 
more and more wide leading to faster escaping orbits, which 
means that the escape period decreases rapidly. We found, 
that the longest escape rates correspond to initial conditions 
near the boundaries between the escape basins and near the 
vicinity of stability islands. On the other hand, the shortest 
escape periods have been measured for the regions without 
sensitive dependence on the initial conditions (basins of es¬ 
cape), that is, those far away from the fractal basin bound¬ 
aries. 

We continue our exploration of the escape process in the 
phase (x,x) plane. The structure of the (x,x) phase plane for 
several values of the energy is shown in Fig. [5] We observe 
a similar behavior to that discussed for the physical (x,y) 
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Fig. 3 Evolution of the percentages of trapped and escaping orbits when varying the energy h (a-left): on the physical (x,y) plane and (b-right): 
on the phase (x,i) plane. 


plane in Fig. [2] The outermost black solid line is the Zero 
Velocity Curve (limiting curve) which is defined as 

f{x,x) = ^i 1 + V(x,y = 0)=h. (11) 

It is worth noticing, that in the phase plane the limiting curve 
is closed but this does not mean that there is no escape. Re¬ 
member, that we decided to choose such perturbation terms 
that produce the escape channels on the physical (x,y) plane 
which is a subspace of the entire four-dimensional (x,y,i,y) 
space of the system. Here we must point out, that this (x,x) 
phase plane is not a Poincare Surface of Section (PSS), sim¬ 
ply because escaping orbits in general, do not intersect the 
y = 0 axis after a certain time, thus preventing us from defy¬ 
ing a recurrent time. A classical Poincare surface of section 
exists only if orbits intersect an axis like y = 0 at least once 
within a certain time interval. Nevertheless, in the case of 
escaping orbits we can still define local surfaces of section 
which help us to understand the orbital behavior of the dy¬ 
namical system. 

Again, we can distinguish in the phase plane fractal re¬ 
gions where we cannot predict the particular escape channel 
and regions occupied by escape basins. These basins are ei¬ 
ther broad well-defined regions, or elongated bands of com¬ 
plicated structure spiralling around the center. We see that 
again for values of energy close to the escape energy there 
is a considerable amount of trapped orbits and the degree 
of fractalization of the phase plane is high. As we proceed 
to higher energy levels however, the rate of trapped orbits 
reduces, the phase plane becomes less and less fractal and 
is occupied by well-defined basins of escape. In Fig. [3J} we 
present the evolution of the percentages of trapped and es¬ 
caping orbits on the phase plane when the value of the en¬ 
ergy h varies. It is observed, that the pattern and the evolu¬ 
tion of the percentages is almost identical to that discussed 


in Fig. [3^ regarding the physical plane. In particular, for 
h = 0.13 about half of the phase plane is covered by initial 
conditions corresponding to trapped orbits, while the escap¬ 
ing orbits share the rest half of the (x,x) plane. At the highest 
energy level studied (h = 0.5), about 60% of the total orbits 
escape through channel 1 and 40% through channel 2; the 
percentage of trapped orbits has already reached the zero 
value from h > 0.3. 

The distribution of the escape times t QSC of orbits on the 
(x,x) plane is shown in Fig. [6] It is evident, that orbits with 
initial conditions inside the exit basins escape from the sys¬ 
tem very quickly, or in other words, they possess extremely 
small escape periods. On the contrary, orbits with initial con¬ 
ditions located in the fractal parts of the phase plane need 
considerable amount of time in order to escape. Another in¬ 
teresting way of measuring the escape rate of an orbit is by 
counting how many intersection the orbit has with the axis 
y = 0 before it escapes. The regions in Fig. [7] are colored 
according to the number of intersections with the axis y = 0 
upwards (y > 0). We observe, that orbits with initial condi¬ 
tions inside the two red basins escape directly without any 
intersection with the y = 0 axis. Furthermore, as the value 
of the energy increases, these red regions grow in relative 
size (proportion of the total area on the phase plane) and for 
high enough energy levels they occupy around 90% of the 
grid. We should also note, that orbits with initial conditions 
located at the vicinity of the stability islands perform numer¬ 
ous intersections with the y = 0 axis before they eventually 
escape to infinity. On the other hand, orbits with initial con¬ 
ditions in the elongated spiral bands need only a couple of 
intersection until they escape. 

The grids in physical (x,y) as well as the phase (x,x) 
plane provide information on the phase space mixing for 
only a fixed value of energy. Henon however, back in the 
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Fig. 4 Distribution of the escape times t esc of the orbits on the (x,y) plane. The darker the color, the larger the escape time. Trapped orbits are 
indicated by gray color. 


60s, considered a plane which provides information about 
regions of stability and regions of escaping orbits using the 
section y = x = 0, y > 0, i.e., the test particle starts on the 
x-axis, parallel to the y-axis and in the positive y-direction. 
Thus, in contrast to the previously discussed grids, only or¬ 
bits with pericenters on the x-axis are included and there¬ 
fore, the value of the energy h is used as an ordinate. Fig. [8] 
shows the structure of the (x, h )-plane when he (0.125,0.5]. 
The boundaries between bounded and unbounded motion 
are now seen to be more jagged than shown in the previous 
grids. In addition, we found in the blow-ups of the diagram 


many tiny islands of stability^] We see, that for low values of 
the energy close to the escape energy, there is a considerable 
amount of trapped orbits inside stability regions surrounding 
by a highly fractal structure. This pattern however changes 
for larger energy levels, where there are no trapped orbits 
and the vast majority of the grid is covered by well-formed 
basins of escape, while fractal structure is confined only near 
the boundaries of the escape basins. 

It is of particular interest to conduct a statistical analysis 
of the escape process in the case of the (x,x) phase plane. 

3 From chaos theory we expect an infinite number of islands of (sta¬ 
ble) quasi-periodic (or small scale chaotic) motion. 
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Fig. 5 The structure of the phase (x,i) plane for several values of the energy h, distinguishing between different escape channels. The color code 
is as follows: Trapped (gray); escape through channel 1 (green); escape through channel 2 (red). 


For this purpose, we shall follow the numerical approach 
used recently in l25l . Our results are shown in Fig. [9ja-d) 
where curve fit approximation versus results from numeri¬ 
cal integration is presented in (a-b) panels. To begin with, 
Fig. shows the proportion of escaping orbits N e /No as a 
function of the energy h. For values of energy beyond the 
escape energy, the majority of orbits escape from the sys¬ 
tem. Our numerical calculation verify that the evolution of 
the proportion of escaping orbits can be approximated by 
the formula 

N e /N 0 (h) = 0.5 [1+tanh (30/1 — 4)], (12) 

proposed in lf25l . In Fig. [9]) we present the evolution of the 
direct escaping orbits Nje/No (by the term direct escaping 


orbits we refer to orbits that escape to infinity immediately 
without any intersection with the y = 0 axis) as a function 
of the energy h. We see, that the amount of direct escaping 
orbits grows rapidly with increasing h and for high energy 
levels (h > 0.5) they take over almost all the phase plane 
(more than 90%). The proportion of direct escapes can be 
given by the approximate formula 

N de /No{h) = —1.7 + 19.24/z — 49.15/z 2 + 42.16/z 3 . (13) 

Moreover, Fig. |9j: depicts the logarithm of the proportion 
of escaping orbits dN n /No, where dN n corresponds to the 
number of escaping orbits after the nth intersection with the 
y = 0 axis upwards (y > 0). It is seen, that the escape time of 
orbits decreases with increasing n. In particular, the escape 
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Fig. 6 Distribution of the escape times t esc of the orbits on the (x,x) plane. The darker the color, the larger the escape time. Trapped orbits are 
indicated by gray color. 


rates are high for relatively small n , while they drop rapidly 
for larger n. Last but not least, we computed the probability 
of escape as a function of the number of intersections for 
various values of the energy. Specifically, the probability is 
defined as 


where N n is the number of orbits that have not yet escaped 
before the nth intersection. The evolution of p n as a func¬ 
tion of n for various energy levels is given in Fig. [9ji. Here 
we have to stress out, that the properties of the probability 
of escape in this system and in other similar systems have 


been studied in detail in previous papers (e.g., es Ha-Ha. 
Furthermore, our numerical calculations regarding the sta¬ 
tistical analysis of the dynamical system in the case where 
two escape channels are present, have found to coincide with 
the corresponding results given in l25l . 


4.2 Case II: Three channels of escape 

We continue our exploration of escapes in a Hamiltonian 
system with three exit channels and escape energy equal to 
1/6. In order to obtain this number of exits in the limiting 
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Fig. 7 Color scale of the escape regions as a function of the number of intersections with the y = 0 axis upwards (y > 0). The color code is as 
follows: 0 intersections (red); 1 intersection (blue); 2 intersections (magenta); 3-10 intersections (orange); >10 intersections (green). The gray 
regions represent stability islands of trapped orbits. 


curve in the (x,y) plane, the perturbation term should be 
V\{x,y) = — v(x 2 /3 —y 2 ) and the corresponding Hamilto¬ 
nian reads 

h 2 = \ b 2 +y 2 +x 2 +y 2 ) -x(x 2 /3-y 2 ) =h. (15) 

#2 manifests a 2^/3 rotation symmetry, but for £ this dis¬ 
crete symmetry is broken. Here we should like to note, that 
the particular type of the perturbation is very similar to that 
of the classical Henon-Heiles Hamiltonian system l36l (in 
fact we changed the position of the v and y variables). We 
made this choice mainly for two reasons: (i) the standard 
Henon-Heiles dynamical system has been studied extensively 
and thoroughly in numerous papers over the last years (e.g., 


m - m isi nu i26i 09)) so, we preferred to work on some¬ 
thing rather different and (ii) in all cases we wanted the 
(x,x) phase plan^J It should be pointed out however, that 
this change in the variables affects only the symmetry, while 
all the measured quantities remain the same as in the clas¬ 
sical Henon-Heiles system. In Fig. [T0] i we see the equipo- 
tential curves of the potential ^ for various values of the 
energy h, while the equipotential corresponding to the en¬ 
ergy of escape h esc is plotted with red color in the same plot. 
Furthermore, the open ZVC at the physical (x,y) plane when 
h = 0.2 > h esc is presented with green color in Fig.fTOfr and 

4 The (jc,jc) phase plane is constructible only if the potential has 
terms with even powers regarding the y variable. 
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Fig. 9 (a-upper left): Evolution of the proportion of escaping orbits N e /N q as a function of the energy h, (b-upper right): Evolution of the proportion 
of directly escaping orbits N^/Nq as a function of the energy h, (c-lower left): Evolution of the logarithmic proportion dN n /No as a function of 
the number of the intersections n, for various values of the energy and (d-lower right): Evolution of the probability p n of escapes as a function of 
n for several energy levels. 



Fig. 10 (a): Equipotential curves of the potential for various values of the energy h when V\ (x,y) = —x(x 2 /3 — y 2 ). The equipotential curve 
corresponding to the energy of escape is shown witn red color; (b): The open ZVC at the physical (x,y) plane when h = 0.2. L\, L 2 and L 3 indicate 
the three unstable Lyapunov orbits plotted in red; (c): Three escaping orbits when h = 0.2. The orbit which escapes from channel 1 is potted with 
green color, the orbit escaping from channel 2 with red color, while blue color is used for the orbits which escapes through channel 3. 
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Fig. 17 Orbital structure of the {x,h )~plane when three channels of es¬ 
cape are present. This diagram gives a detailed analysis of the evolution 
of the trapped and escaping orbits of the dynamical system when the 
parameter h changes. The color code is as in Fig. m 


the three channels of escape are shown. In the same figure, 
the three unstable Lyapunov orbits L \, L 2 and L 3 are denoted 
using red color. Fig. [TO]: depicts with different colors three 
orbits, one escaping from channel 1 , one from channel 2 and 
the other from channel 3, when h = 0.2. 


In this case, we shall investigate the escape properties of 
unbounded motion of test particles for values of energy in 
the set h = {0.17,0.18,0.19,0.20,0.22,0.24,0.26,0.28,0.30}. 
We begin with initial conditions of orbits in the physical 
(x,y) plane. The orbital structure of the physical plane for 
different values of the energy h is show in Fig. 11 Again, 
following the approach of the previous case, each initial con¬ 
dition is colored according to the escape channel through 
which the particular orbit escapes. Stability islands on the 
other hand, filled with initial conditions of orbits which do 
not escape are indicated as gray regions. We observe, that 
things are quite similar to that discussed previously in Fig. 
[2] In fact, for energy levels very close to the escape energy, 
the central region of the plot is highly fractal and it is also oc¬ 
cupied by several stability islands. However, as we increase 
the value of the energy the regions of regular trapped orbits 
are reduced, the physical plane becomes less and less fractal 
and well-defined basins of escape emerge. 


The evolution of the percentages of trapped and escap¬ 
ing orbits on the physical (x, y) plane when the value of the 
energy h varies is presented in Fig.[T2]i. It is seen, that when 
h = 0.17, that is the first investigated energy level above 
the escape energy, escaping orbits through channels 2 and 3 
share the same percentage (around 32%), escapers through 


channel 1 have a slightly elevated percentage (around 34%), 
while trapped orbits possess a very low rate corresponding 
only to 4% of the physical plane. Once more, as we increase 
the value of the energy the rate of trapped orbits decreases 
and eventually vanishes for h > 0.2. Furthermore, we ob¬ 
serve that the percentage of escaping orbits through chan¬ 
nel 2 grows with increasing energy and for h > 0.4 it seems 
to saturate around 44%. The percentage of escaping orbits 
through channel 3 on the other hand, exhibits a slow but 
constant decrease, while the rate of escaping orbits through 
exit 1 after small fluctuations it saturates around 32% for 
h > 0.4. In general terms, we may conclude that throughout 
the energy range studied, the majority of orbits in the phys¬ 
ical (x,y) plane choose to escape through channel 2 , while 
exit 3 seems to be the least favorable among the escape chan¬ 
nels. 


The following Fig. [13] shows how the escape times t esc of 
orbits are distributed on the (x,y) plane. Light reddish colors 
correspond to fast escaping orbits, dark blue/purpe colors in¬ 
dicate large escape periods, while gray color denote trapped 
orbits. This grid representation of the physical plane gives 
us a much more clearer view of the orbital structure and es¬ 
pecially about the trapped orbits. In particular, we see that 
for h = 0.2 we have the last indication of stability islands, 
as for all higher energy levels studied all orbits escape, thus 
defying basins of escape. 

Our exploration continuous in the phase (x, x) plane. The 
structure of the (x,x) phase plane for different values of the 
energy is shown in Fig. [14] We observe a similar behavior to 
that discussed for the physical (x,y) plane in Fig. 11 Again, 
we can distinguish in the phase plane fractal regions where 
the prediction of the particular escape channel is impossible 
and regions occupied by escape basins. It is interesting to 
note, that the limiting curve (ZVC) is open at the right part 
due to the x 3 term entering the perturbation function. The 
rich fractal structure of the phase space shown in the grids of 
Fig. [14] implies that our system has also a strong topological 
property, which is known as the Wada property. This special 
topological property has been identified and studied in sev¬ 
eral dynamical systems (e.g., a sni mi) and it is a typical 
property in open Hamiltonian systems with three or more 
escape channels. An escape basin is a Wada basin if any 
boundary point also belongs to the boundary of at least two 
other basins |QT| 40|. It is seen in Fig. [14] that for h > 0.25 
all the KAM regime vanishes 0 and therefore, all the initial 
conditions of orbits escape through one of the exits. 


It is evident from Fig.fTZfr where the evolution of the per¬ 
centages of trapped and escaping orbits on the phase plane 
as a function of the value of the energy h is presented, that 
the pattern has many differences comparing to that discussed 
previously in Fig. [I2fc ; only the percentage of trapped or¬ 
bits exhibits similar behavior. To begin with, we observe 
that for h = 0.17 more than half of the phase plane (around 
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Fig. 11 The structure of the physical (jt,y) plane for several values of the energy h, distinguishing between different escape channels. The color 
code is as follows: Trapped (gray); escape through channel 1 (green); escape through channel 2 (red); escape through channel 3 (blue). 


55%) corresponds to initial conditions of orbits that escape 
through channel 1, while orbits escaping through exits 2 and 
3 share about 44% of the grid. As the value of the energy in¬ 
creases and we move away from the escape energy it is seen, 
that the rate of orbits escaping through exit 1 increases and 
always dominates, while on the hand, the percentages of or¬ 
bits escaping through channels 2 and 3 drop. At the highest 
energy level studied (h = 0.5), about 70% of the total or¬ 
bits escape through channel 1, about 20% through channel 2 
and only 10% through channel 3. Thus, one may reasonably 
conclude that throughout the energy range studied, the vast 
majority of orbits in the phase (x,i) plane choose to escape 
through channel 1, while channels 2 and 3 are much less 
likely to be chosen. 


Fig. |T5] shows the distribution of the escape times t Qsc 
of orbits on the (v,i) plane. It is evident, that orbits with 
initial conditions inside the exit basins escape from the sys¬ 
tem after short time intervals, or in other words, they pos¬ 
sess extremely small escape periods. On the contrary, or¬ 
bits with initial conditions located in the fractal parts of the 
phase plane need considerable amount of time in order to 
find one of the exits and escape. We see, that for h > 0.2 
there is no indication of stability islands corresponding to 
trapped orbits. In another point of view, Fig. |T6] shows the 
regions of the phase plane which are now colored accord¬ 
ing to the number of intersections the orbits perform with 
the axis y = 0 upwards (y > 0). The red regions denote ini¬ 
tial conditions of orbits that escape directly from the system 
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Fig. 12 Evolution of the percentages of trapped and escaping orbits when varying the energy h (a-left): on the physical (x,y) plane and (b-right): 
on the phase (x,x) plane. 


without ever intersecting the y = 0 axis. The proportion of 
the total area on the phase plane occupied by these regions 
of direct escapes grows with increasing energy and for high 
enough energy levels they occupy more than 70% of the 
grid. In Fig. [IT] we present the structure of the (x,h )-plane 
when he ( 1/6,1/2]. It is seen, that trapped orbits exist only 
at low energies very close to the escape energy (h < 0 . 22 ), 
while for larger energy levels all the orbits escape to infinity. 
Once more, highly fractal structure is observed near the sta¬ 
bility islands of regular motion, while the degree of fractal- 
ization, or in other words the unpredictability of the system, 
reduces significantly where there are no trapped orbits and 
well-defined basins of escape cover the vast majority of the 
(x,h) -plane. 

At this point, we shall follow the approach discussed in 
subsection 4.2 in order to perform a statistical analysis of the 
escape process in the case of the (x,x) phase plane for the 
Hamiltonian system with three channels of escape. Fig. |~j~ 8 } i 
shows the proportion of escaping orbits N e /No as a function 
of the energy h. For values of energy beyond the escape en¬ 
ergy, more than 95% of the total orbits escape from the sys¬ 
tem. According to our numerical calculation, the evolution 
of the proportion of escaping orbits can be approximated by 
the formula 


N e /N 0 (h) = 0.5 [1+tanh ( 49/1 — 6.5)]. (16) 

In Fig. [TSfr we present the evolution of the direct escaping 
orbits Nde/No as a function of the energy h. We see, that 
the amount of direct escaping orbits grows rapidly with in¬ 
creasing h and for high energy levels (h > 0.5) they populate 
about 80% of the phase plane. The proportion of direct es¬ 
capes can be given by the approximate formula 

N de /N 0 {h) = —2.2 + 22.86/j — 59.24/r 2 + 51 ,42h 2 . (17) 


Furthermore, Fig. [18p depicts the logarithm of the propor¬ 
tion of escaping orbits dN n /No , as a function of the inter¬ 
sections with the y = 0 axis upwards an orbit performs be¬ 
fore it escapes. We observe, that the escape time of orbits 
decreases with increasing n. In particular, the escape rates 
are high for relatively small n , while they drop rapidly for 
larger n. Finally, we calculated the probability of escape as a 
function of the number of intersections for various values of 
the energy. The evolution of p n as a function of n for various 
energy levels is shown in Fig.[l 8 H. 


4.3 Case III: Four channels of escape 

The last case under investigation is a Hamiltonian system 
with four channels of escape. In order to obtain this number 
of exits in the limiting curve in the physical (x, y) plane, we 
chose the perturbation term V\ (x,y) = — x 2 y 2 and the corre¬ 
sponding Hamiltonian is he following 

H 3 = 1 (x 2 +y 2 +x 2 +y 2 )-x 2 y 2 = h. (18) 

The Hamiltonian H 3 is invariant under x — x and/or y —>> 
—y. The escape mechanism in this particular Hamiltonian 
system with the four escape channels and escape energy 
equals to 1/4 has already been examined (e.g., I20ll23ll24l 
i38'.B21 ). In Fig.[T 9 }i, we present the equipotential curves of 
the potential for various values of the energy h, while 
the equipotential corresponding to the energy of escape h esc 
is plotted with red color in the same plot. In addition, the 
open ZVC at the physical (x,y) plane when h = 0.3 > h esc is 
given with green color in Fig. [T9fr , while the four channels 
of escape are also shown. In the same figure, the four unsta¬ 
ble Lyapunov orbits L\, L 2 , Z /3 and L 4 are denoted using red 
color. In Fig.|T9jc we plotted with different colors four orbits, 
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Fig. 13 Distribution of the escape times t Q sc of the orbits on the (x,y) plane. The darker the color, the larger the escape time. Trapped orbits are 
indicated by gray color. 


one escaping from channel 1, one from channel 2, one from 
channel 3 and the last one from channel 4, when h = 0.3. 


The escape properties and mechanism of unbounded mo¬ 
tion of test particles for values of energy in the set h = 
{0.26,0.28,0.30,0.35,0.40,0.45,0.50,0.55,0.60} will be ex¬ 
amined. We begin, as usual, with initial conditions of or¬ 
bits in the physical (v,y) plane. Fig. 20 shows the orbital 
structure of the physical plane for different values of the en¬ 
ergy h. Again, following the same approach of the previous 
cases, each initial condition is colored according to the es¬ 
cape channel through which the particular orbit escapes. Ar¬ 
eas corresponding to trapped orbits on the other hand, are 


indicated as gray regions. It is evident, that the structure of 
the (x,y) plane differs significantly with respect to the plots 
shown previously in Figs.[2]and[TT] We see, that for values of 
energy very close to the escape energy almost all the central 
region of the grid is covered by initial conditions of trapped 
orbits, while escaping orbits exist only near the four exits. 
However, with increasing energy the area on the physical 
plane occupied by trapped orbits reduces and several basins 
of escape begin to emerge. At the highest energy level stud¬ 
ied (h = 0.6), there is no indication of trapped motion and all 
orbits escape to infinity through one of the four escape chan¬ 
nels. We also observe, the existence of well-formed basins 
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Fig. 14 The structure of the phase (x.x) plane for several values of the energy h, distinguishing between different escape channels. The color code 
is as follows: Trapped (gray); escape through channel 1 (green); escape through channel 2 (red); escape through channel 3 (blue). 


of escape, while the central region of the grid still remains 
highly fractal. Here we should like to note, that in general 
terms, throughout the energy range the structure of the phys¬ 
ical plane (x,y) is symmetrical with respect to the x = 0 axis. 

It is of particular interest to monitor the evolution of the 
percentages of trapped and escaping orbits on the physical 
(x,y) plane when the value of the energy h varies. A diagram 
depicting this evolution is presented in Fig. |2lfc. We see, that 
for h = 0.26, that is an energy level just above the escape 
energy, about 50% of the physical plane is covered by ini¬ 
tial conditions of trapped orbits. As the value of the energy 
increases however, the rate of trapped orbits drops rapidly 
and eventually at h = 0.6 it vanishes. We also observe, that 
the evolution of the percentages of orbits escaping through 


channels 1 and 3 coincide with the evolution of the percent¬ 
ages escaping through channels 2 and 4, respectively. We 
anticipated this behaviour of the escape percentages, which 
is a natural result of the symmetrical structure of the (x,y) 
plane. It is seen, that initially (h = 0.26) all rates of escaping 
orbits coincide at about 14%. Then, with increasing energy 
the rates of escaping orbits increase and also start to diverge. 
At the highest energy studied, escaping orbits through chan¬ 
nels 1 and 2 share about 65% of the physical plane, while 
escaping orbits trough channels 3 and 4 occupy the remain¬ 
ing 35% of the grid. Therefore, one may reasonably con¬ 
clude that in general terms, throughout the range of the val¬ 
ues of the energy studied, the majority of orbits in the physi- 
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Fig. 15 Distribution of the escape times t esc of the orbits on the (x.x) plane. The darker the color, the larger the escape time. Trapped orbits are 
indicated by gray color. 


cal (x,y) plane choose to escape either through channel 1 or 
channel 2. 

The distribution of the escape times t esc of orbits on the 
physical plane is given in Fig. [22} Light reddish colors cor¬ 
respond to fast escaping orbits, dark blue/purpe colors indi¬ 
cate large escape periods, while gray color denote trapped 
orbits. Here, we have a better view regarding the amount of 
trapped orbits. Indeed, we see that for h = 0.6 all orbits es¬ 
cape from the system. Moreover we observe, that orbits with 
initial conditions close to the area occupied by trapped orbits 
have significantly large escape periods, while on the other 


hand, orbits located near the escape channels escape very 
quickly having escaping rates of about two orders smaller. 


We continue our investigation to the phase (x, x) plane, 
the structure of which for different values of the energy is 
presented in Fig. 23 One may observe, that for h < 0.3 most 
of the phase plane is covered by a vast region corresponding 
to trapped orbits, while only two small islands of initial con¬ 
ditions of escaping orbits exit. However, as the value of the 
energy increases and we move far away for the escape en¬ 
ergy, the extent of these two islands grows and for h > 0.35 
the trapped orbits are mainly confined to the central region 
of the phase plane. At the same time, small elongated spiral 




















Escapes in Hamiltonian systems with multiple exit channels: Part I 


19 




-0.5 0.0 0.5 1.0 1.5 

h = 0.20 




0.5 

x 

-0.5 



h = 0.22 


1.5 

1.0 

0.5 

-0.5 

- 1.0 

-1.5 


x 



1.5 

1.0 

0.5 

-0.5 

- 1.0 

-1.5 




1.5 

1.0 

0.5 

-0.5 

- 1.0 

-1.5, 


1.5 

1.0 

0.5 

-0.5 

- 1.0 

-1.5 


Fig. 16 Color scale of the escape regions as a function of the number of intersections with the y = 0 axis upwards (y > 0). The color code is as 
follows: 0 intersections (red); 1 intersection (blue); 2 intersections (magenta); 3-10 intersections (orange); >10 intersections (green). The gray 
regions represent stability islands of trapped orbits. 


basins of escape emerge inside the fractal region which sur¬ 
rounds the area of trapped orbits. Furthermore, at very high 
energy levels (h > 0.55) we see that trapped orbits disap¬ 
pear completely from the grid and the two main basins of 
escape take over the vast majority of the phase plane, while 
the elongated escape basins remain confined to the central 
region. As we noticed previously when discussing the phys¬ 
ical (v,y) plane, there is also a symmetry in the phase plane. 
In particular, throughout the energy range the structure of 
the phase plane (v, x) is somehow symmetrical (not with the 
strick sense) with respect to the x = 0 axis. 

The evolution of the percentages of trapped and escaping 
orbits on the phase plane as a function of the value of the en¬ 


ergy h is given in Fig.|2ljx For h = 0.26, we see that trapped 
orbits dominate the phase plane as they occupy about 90% 
of the gird. However as usual, with increasing energy the 
dominance of trapped orbits deteriorates rapidly due to the 
increase of the rates of escaping orbits which form basins of 
escape. We observe, that once more as in Fig.[2T[r the evolu¬ 
tion of the percentages of orbits escaping through channels 1 
and 3 coincides with the evolution of the percentages escap¬ 
ing trough channels 2 and 4, respectively. The percentages of 
all types of escaping orbits increase but with different rates 
and for h > 0.35 they overwhelm the amount of trapped or¬ 
bits. In particular, we see that the percentages of orbits es¬ 
caping through exits 1 and 2 are always higher than those 
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Fig. 18 (a-upper left): Evolution of the proportion of escaping orbits N e /No as a function of the energy h, (b-upper right): Evolution of the 
proportion of directly escaping orbits N^/Nq as a function of the energy h, (c-lower left): Evolution of the logarithmic proportion dN n /No as a 
function of the number of the intersections n, for various values of the energy and (d-lower right): Evolution of the probability p n of escapes as a 
function of n for several energy levels. 



Fig. 19 (a): Equipotential curves of the potential |2J for various values of the energy h when Vi(x,y) = —x 2 y 2 . The equipotential curve corre¬ 
sponding to the energy of escape is shown with red color; (b): The open ZVC at the physical (x,y) plane when h = 0.3. Li, L 2 , L 3 and L 4 indicate 
the four unstable Lyapunov orbits plotted in red; (c): Four escaping orbits when h — 0.3. The orbit which escapes from channel 1 is potted with 
green color, the orbit escaping from channel 2 with red color, the one from channel 3 with blue, while orange color is used for the orbits which 
escapes through channel 4. 
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Fig. 20 The structure of the physical (%,y) plane for several values of the energy h, distinguishing between different escape channels. The color 
code is as follows: Trapped (gray); escape through channel 1 (green); escape through channel 2 (red); escape through channel 3 (blue); escape 
through channel 4 (orange). 


corresponding to orbits escaping through channels 3 and 4. 
Moreover, the rates of exits 1 and 2 increase constantly and 
at the highest energy level studied (h = 0.6) the share about 
90% of the entire phase plane. On the other hand, the per¬ 
centages of exits 3 and 4, even though they also grow with 
increasing energy, they always possess significantly smaller 
values than exits 1 and 2 and for h > 0.4 they seem to sat¬ 
urate around 5%. Thus, we may conclude that the vast ma¬ 
jority of orbits in the phase (x,i) plane exhibit clear sings 
of preference through exits 1 and 2, while channels 3 and 4 
have considerable less probability to be chosen. 

Taking into account that for low values of the energy 
there is a considerable amount of trapped motion in the phase 


plane, we decided to use the SALI method in order to dis¬ 
tinguish between regular and chaotic trapped orbits. In Fig. 
[24} i we present the phase plane for h = 0.26, where each 
initial condition is plotted according to the regular (cyan) or 
chaotic (magenta) character of the orbit, while white areas 
correspond to escaping orbits. It is seen, that the vast ma¬ 
jority of the trapped orbits are regular however, a thin layer 
composed of chaotic trapped orbits is also present. There¬ 
fore, a natural and very important question arises: do these 
chaotic bounded orbits remain trapped forever? Remember, 
that in the current investigation we set the maximum time of 
the numerical integration to be equal to 10 5 time units. We 
suspect, that all these trapped chaotic orbits will eventually 
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Fig. 21 Evolution of the percentages of trapped and escaping orbits when varying the energy h (a-left): on the physical (x,y) plane and (b-right): 
on the phase (x,x) plane. 



Fig. 27 Orbital structure of the (x, h) -plane when four escape channels 
are present. This diagram gives a detailed analysis of the evolution of 
the trapped and escaping orbits of the dynamical system when the pa¬ 
rameter h changes. The color code is as in Fig.|20| 


escape from the system if they have enough time to evolve. 
Thus, in order to shed some light to this issue, we let the time 
running and we integrated these orbits until they escape. Our 
numerical calculations revealed, that these orbits are in fact 
super sticky orbits which possess extremely high escape pe¬ 
riods up to 3.5 x 10 6 time units. A characteristic example 
of such a super sticky orbit is given in Fig. [24fr , where we 
monitor the time evolution of R 2 = x 2 -fy 2 . The horizon¬ 
tal, red, dashed line at 1.26 approximates the position of the 
unstable Lyapunov orbits at the four exits, while the ver¬ 


tical, blue, dashed line denotes the initial integration time 
(10 5 time units). We see, that the particular orbit escapes 
through channel 2 after a time interval of about 851000 time 
units which is more than 8.5 times the initial integration pe¬ 
riod. Our additional computations indicate, that these super 
sticky orbits correspond to less than 10% of the total (regular 
plus chaotic) trapped orbits so, using 10 5 time units for the 
numerical integration and counting them as trapped, even 
though they escape after vast time intervals, does not have a 
huge impact in our results. 


The following Fig. 25 shows the distribution of the es¬ 
cape times t Qsc of orbits on the (x,i) plane. It is clear, that 
orbits with initial conditions inside the exit basins escape 
to infinity after short time intervals, or in other words, they 
possess extremely small escape periods. On the contrary, or¬ 
bits with initial conditions located in the fractal parts of the 
phase plane need considerable amount of time in order to 
find one of the four exits and escape. It is seen, that at the 
highest energy level studied (h = 0.6) there is no indication 
of bounded motion and all orbits escape to infinity sooner or 
later. 


In Fig. [26} we reconstructed the grids in the phase plane 
using different color code and now the regions of the phase 
plane are colored according to the number of intersections 
the orbits perform with the axis y = 0 upwards (y > 0). 
Specifically, red regions correspond to initial conditions of 
orbits that escape directly from the system without any in¬ 
tersection with the y = 0 axis. We see, that the proportion of 
the total area on the phase plane occupied by orbits which 
escape directly from the system grows rapidly with increas¬ 
ing energy and for h > 0.6 they occupy more than 90% of 
the entire grid. Fig. 27 depicts the structure of the (x,h)- 
plane when h E (0.25,1]. At low energy levels, one may ob¬ 
serve, three important issues: (i) the vast majority of oribis 
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Fig. 22 Distribution of the escape times t esc of the orbits on the (%,y) plane. The darker the color, the larger the escape time. Trapped orbits are 
indicated by gray color. 


are trapped, (ii) the structure of the (x,h )-plane exhibits a 
high degree of fractalization and (iii) basins of escape are 
present only at the outer parts of the grid. However, when 
h > 0.6, that is when trapped orbits cease to exist we see 
that the fractal structure disappears and all the (x,h }-plane 
is coved by well-defined basins of escape. It should also be 
pointed out, that the structure of the (x,h )-plane is symmet¬ 
rical with respect to the v = 0 axis. 

Before closing this section, we would like to perform 
a statistical analysis of the escape process in the case of 
the (x,x) phase plane for the Hamiltonian system with four 
channels of escape. The proportion of escaping orbits N e /No 


as a function of the energy h is presented in Fig.|28p. We see 
that for h > 0.45, more than 90% of the total orbits escape 
from the system. Our numerical computations suggest, that 
the evolution of the proportion of escaping orbits can be ap¬ 
proximated by the formula 


N e /N 0 {h ) = 0.5 [1 +tanh(15.85/* -4.93)]. (19) 

Furthermore, in Fig. [28]3 we present the evolution of the di¬ 
rect escaping orbits N^ e /No as a function of the energy h. As 
it was found in the previously examined cases, the amount 
of direct escaping orbits grows rapidly with increasing h and 
for high energy levels (h > 0.6) they populate more than 
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Fig. 23 The structure of the phase (v,x) plane for several values of the energy h, distinguishing between different escape channels. The color code 
is as follows: Trapped (gray); escape through channel 1 (green); escape through channel 2 (red); escape through channel 3 (blue); escape through 
channel 4 (orange). 


90% of the phase plane. The proportion of direct escapes 
can be given by the approximate formula 


N de /N 0 {h) = —1.496 + 7.585/i — 5.938/i 2 


( 20 ) 


The evolution of the logarithm of the proportion of escap¬ 
ing orbits dN n /No , as a function of the intersections with 
the y = 0 axis upwards an orbit performs before it escapes is 


given in Fig.|28p. One may observe, that the escape time of 
orbits decreases with increasing n. Being more precise, the 
escape rates are high enough for relatively small number of 
intersections n , while they fall rapidly for larger n. Finally, 
we computed the probability of escape as a function of the 
number of intersections for various values of the energy h. 


Our results are shown in Fig.|28]i, where we present the evo¬ 
lution of p n as a function of n for various energy levels. 


5 Conclusions and discussion 

The main objective of this work was to review but also nu¬ 
merically investigate even further the escape properties of 
orbits in a dynamical system of two-dimensional coupled 
perturbed harmonic oscillators, which is a characteristic ex¬ 
ample of open Hamiltonian systems. The key feature of this 
type of Hamiltonians is that they have a finite energy of 
escape. In particular, for energies smaller than the escape 
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Fig. 24 (a-left): The structure of the phase (x,x) plane for h = 0.26, distinguishing between trapped regular orbits (cyan), trapped chaotic orbits 
(magenta) and escaping orbits (white), (b-left): Time evolution of R 2 = x 2 +y 2 for a super sticky orbit when h = 0.26. The horizontal, red, dashed 
line at 1.26 approximates the position of the unstable Lyapunov orbits at the four exits, while the vertical, blue, dashed line denotes the initial 
integration interval of 10 5 time units. 


value, the equipotential surfaces are close and therefore es¬ 
cape is impossible. For energy levels larger than the escape 
energy however, the equipotential surfaces open and several 
channels of escape appear through which the particles can 
escape to infinity. Here we should emphasize, that if a test 
particle has energy larger than the escape value, this does not 
necessarily mean that the particle will certainly escape from 
the system and even if escape does occur, the time required 
for an orbit to cross a Lyapunov orbit and hence escape to in¬ 
finity may be vary long compared with the natural crossing 
time. The function containing the perturbation terms affects 
significantly the structure of the equipotential curves and de¬ 
termines the exact number of the escape channels. We chose 
such forms of perturbations and divided our study into three 
cases with respect to the number of the escape channels. 

Since a distribution function of the system was not avail¬ 
able so as to use it for extracting the different samples of 
orbits, we had to follow an alternative path. We defined for 
each set of values of the energy, dense grids of initial con¬ 
ditions regularly distributed in the area allowed by the value 
of the energy in both the physical and the phase space. In 
both cases, the density of the grids was controlled in such 
a way that always there are about 50000 orbits to be exam¬ 
ined. For the numerical integration of the orbits in each grid, 
we needed roughly between 1 minute and 3 days of CPU 
time on a Pentium Dual-Core 2.2 GHz PC, depending both 
on the amount of trapped orbits and on the escape rates of 
orbits in each case. For each initial condition, the maximum 
time of the numerical integration was set to be equal to 10 5 
time units however, when a particle escapes the numerical 
integration is effectively ended and proceeds to the next ini¬ 
tial condition. 


The structure of both the physical (x,y) and phase (x,x) 
space has been explored for several values of the energy h. 
We managed to distinguish between trapped (non-escaping) 
and escaping orbits and we also located the basins of es¬ 
cape leading to different exit channels, finding correlations 
with the corresponding escape times of the orbits. Among 
the escaping orbits, we separated between those escaping 
fast or late from the system. Our extensive numerical cal¬ 
culations strongly suggest, that the overall escape process is 
very dependent on the value of the total orbital energy. We 
also performed a statistical analysis in each case, relating the 
proportion of escaping and directly escaping orbits with the 
value of the energy. In the same vein, the evolution of the 
proportion of escaping orbits and the corresponding proba¬ 
bility, as functions of the nth intersection with the y = 0 axis 
upwards was also presented. 

The main numerical results of our investigation can be 
summarized as follows: 

1. In all three cases studied, areas of trapped orbits and re¬ 
gions of initial conditions leading to escape in a given 
direction (basins of escape), were found to exist in both 
the physical and the phase space. The several escape 
basins are very intricately interwoven and they appear 
either as well-defined broad regions or thin elongated 
spiral bands. Regions of trapped orbits first and foremost 
correspond to stability islands of regular orbits where a 
third integral of motion is present. 

2. A strong correlation between the extent of the basins of 
escape and the value of the energy h was found to exists. 
Indeed, for low values of h the structure of both physical 
and phase space exhibits a large degree of fractalization 
and therefore the majority of orbits escape choosing ran- 
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Fig. 25 Distribution of the escape times t Qsc of the orbits on the (x.x) plane. The darker the color, the larger the escape time. Trapped orbits are 
indicated by gray color. 


domly escape channels. As the value of h increases how¬ 
ever, the structure becomes less and less fractal and sev¬ 
eral basins of escape emerge. The extent of these basins 
of escape is more prominent at high energy levels, where 
they occupy about 90% of the entire area on the grids. 

3. It was found, that for energy levels slightly above the es¬ 
cape energy the majority of the escaping orbits have con¬ 
siderable long escape rates (or escape periods), while as 
we proceed to higher energies the proportion of fast es¬ 
caping orbits increases significantly. This phenomenon 
can be justified, if we take into account that with increas¬ 
ing energy the exit channels on the equipotential curves 


become more and more wide thus the test particles can 
find easily and faster one of the exits and escape to in¬ 
finity. 

4. We observed, that in several exit regions the escape pro¬ 
cess is highly sensitive dependent on the initial condi¬ 
tions, which means that a minor change in the initial con¬ 
ditions of an orbit lead the test particle to escape through 
another exit channel. These regions are the opposite of 
the escape basins, are completely intertwined with re¬ 
spect to each other (fractal structure) and are mainly lo¬ 
cated in the vicinity of stability islands. This sensitiv¬ 
ity towards slight changes in the initial conditions in the 
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Fig. 26 Color scale of the escape regions as a function of the number of intersections with the y = 0 axis upwards (y > 0). The color code is as 
follows: 0 intersections (red); 1 intersection (blue); 2 intersections (magenta); 3-10 intersections (orange); >10 intersections (green). The gray 
regions represent stability islands of trapped orbits. 


fractal regions implies, that it is impossible to predict 
through which exit the particle will escape. 

5. Our calculations revealed, that the escape times of orbits 
are directly linked to the basins of escape. In particular, 
inside the basins of escape as well as relatively away 
from the fractal domains, the shortest escape rates of the 
orbits had been measured. On the other hand, the longest 
escape periods correspond to initial conditions of orbits 
either near the boundaries between the escape basins or 
in the vicinity of stability islands. 

6. In the case where the perturbation term creates four chan¬ 
nels of escape in the physical space, we found that a 
small portion of chaotic orbits with initial conditions 


close to the outermost KAM islands remain trapped in 
the neighbourhood of these islands for vast time inter¬ 
vals having sticky periods which correspond to hundreds 
of thousands time units. On the contrary in systems with 
two and three exit channels all non-escaping orbits are 
regular, while all tested chaotic orbits escape to infinity 
within the predefined integration time. 

We hope that the present review analysis and the corre¬ 
sponding numerical results to be useful in the active field 
of open Hamiltonian systems which may have implications 
in different aspects of chaotic scattering with applications 
in several areas of physics. In Part II of our investigation, 
we shall try to reveal the escape properties of orbits in dy- 
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Fig. 28 (a-upper left): Evolution of the proportion of escaping orbits N e /No as a function of the energy h, (b-upper right): Evolution of the 
proportion of directly escaping orbits N^/Nq as a function of the energy h, (c-lower left): Evolution of the logarithmic proportion dN n /No as a 
function of the number of the intersections n, for various values of the energy and (d-lower right): Evolution of the probability p n of escapes as a 
function of n for several energy levels. 


namical systems with n (n > 5) channels of escape in the 
physical space. Furthermore, it is in our future plans to ex¬ 
pand our exploration in other more complicated potentials, 
focusing our interest in revealing the escape process of or¬ 
bits of stars in realistic galactic systems (i.e., star clusters, 
binary stellar systems, rotating galaxies leaking stars, etc). 
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